function [tides_out] = plot_compare_obsROMS_tides(basename,romsgaugedir,obsdir,gaugeNames,constit)
% plot_compare_obsROMS_tides.m
%
% [tides_out] = plot_compare_obsROMS_tides(basename,romsgaugedir,obsdir,gaugeNames,constit)
%
% A script to make some specific plots of tidal parameters for comparing
% model output to observations.  Various input parameters are be defined
% below and passed to the companion function 'compare_obsROMS_tides.m'.
% Refer to compare_obsROMS_tides.m for a description of the statistical
% comparisons.
%
% Plots of the following are made:
% 1) a map of tide gauge station locations
% 2) labeled bar plots of amplitude and phase (obs & model) for each
%    specified constituent
% 3) a map showing a transect into Puget Sound and Strait of Georgia also
%    with the tide gauge locations 
% 4) line plots of amplitude and phase (obs & model) for each constituent
%    as a function of inland distance relative to Neah Bay, WA (some sites
%    such as Tofino, BC and Bamfield, BC are omitted from this plot)
%
%
% INPUT:  (optional, can be specifed in code below)
%
%     'basename'      -  model run basename
%     'romsgaugedir'  -  path location of ROMS tide gauge extractions
%     'obsdir'        -  path location of observational tide gauge records
%                        currently assumes the following structure:
%                          obsdir/NOAA/NOAA_year
%                          obsdir/Canada/CAN_year
%     'gaugeNames'    -  vector of strings of tide gauge station names
%                        (e.g., 'La Push','Seattle,'Vancouver', etc.)
%     'constit'       -  vector of tidal constituent names for which to
%                        perform the analysis (e.g., 'M2','S2','K1','O1', etc.)
%
%
% R. McCabe 06/2014


%% Set global paths and directories
addpath('../../alpha/');
Tdir = toolstart;
addpath(Tdir.tools);
addpath(Tdir.data);
addpath(genpath([Tdir.tools 'validation_tools']));


%% Define the various input parameters to pass in to the comparison routine
% Set default path locations, basename, gaugeNames, constit if none were provided
%
% Set the run basename
if ~exist('basename','var')
    basename = 'T40final_base_2005_noramp';
end
% Location of ROMS tide gauge extractions
if ~exist('romsgaugedir','var')
    romsgaugedir = '/data1/rmccabe/runs/PostProcess/ROMS_tides/';
end
% Parent directory for the observations
if ~exist('obsdir','var')
    %obsdir = '/data1/rmccabe/tools_data/tide_gauge';
    obsdir = [Tdir.data 'obs_data/tide_gauge'];
end
% Define which tidal constituents to compute
if ~exist('constit','var')
    constit = {'M2';'S2';'K1';'O1'};
end
% Set list of tide gauge locations to extract - comment out the ones you don't want!
if ~exist('gaugeNames','var')
    gaugeNames = {'Charleston'
        'South Beach'
        'Garibaldi'
        'Astoria'
        'Toke Point'
        %'Westport' % exists for 2006 onward
        'La Push'
        'Bamfield'
        'Tofino'
        'Neah Bay'
        'Port Angeles'
        'Victoria'
        'Port Townsend'
        'Friday Harbor'
        'Patricia Bay'
        'Cherry Point'
        'Vancouver'
        'Point Atkinson'
        %'Campbell River' % just outside of Cascadia domain
        'Seattle'
        'Tacoma'};
end

% Fontsize for plotting
fs1 = 12;


%% Perform the tidal analysis and comparisons
[tides_out] = compare_obsROMS_tides(basename,romsgaugedir,obsdir,gaugeNames,constit);


%% Plot a map with the stations
figure; hold on;
set(gcf,'position',[100 290 400 630]);
plot_WAcoast('all','color',[0.5 0.5 0.5]);
axis([-127 -122 43 50.5]);
Z_dar;

for ii = 1:length(tides_out.station_name)
    plot(tides_out.lon(ii),tides_out.lat(ii),'ko','markerfacecolor','g','markersize',6);
    
    % label station names
    switch tides_out.station_name{ii}
        case {'Garibaldi','Toke Point','Westport','La Push','Port Townsend','Vancouver','Campbell River','Seattle','Tacoma'}
            text(tides_out.lon(ii)+0.2,tides_out.lat(ii)-0.15,tides_out.station_name(ii));
        case 'Neah Bay'
            text(tides_out.lon(ii)-1.2,tides_out.lat(ii)+0.1,tides_out.station_name(ii));
        case 'Point Atkinson'
            text(tides_out.lon(ii)+0.2,tides_out.lat(ii)+0.1,tides_out.station_name(ii));
        case 'Friday Harbor'
            text(tides_out.lon(ii)+0.2,tides_out.lat(ii)-0.05,tides_out.station_name(ii));
        case 'Victoria'
            text(tides_out.lon(ii)+0.1,tides_out.lat(ii)-0.07,tides_out.station_name(ii));
        case 'Astoria'
            text(tides_out.lon(ii),tides_out.lat(ii)-0.2,tides_out.station_name(ii));
        case 'Port Angeles'
            text(tides_out.lon(ii)-0.75,tides_out.lat(ii)+0.07,tides_out.station_name(ii));
        otherwise
            text(tides_out.lon(ii)+0.2,tides_out.lat(ii),tides_out.station_name(ii));
    end
end
xlabel('Longitude [deg]','fontsize',12);
ylabel('Latitude [deg]','fontsize',12);
%print(gcf,'-dpng',[outdir 'TideStationLocations'])


%% Plot amplitude and phase bar plots for each tidal constituent
for ii = 1:length(constit)
    constname = char(constit(ii)); % the current constituent
    
    figure; hold on;
    set(gcf,'position',[100 300 900 450]);
    %set(gcf,'position',[101 297 902 468],'paperpositionmode','auto');
    
    subplot_close(2,1,1); hold on;
    bar(1:length(tides_out.lat),tides_out.(constname).amp);
    set(gca,'xlim',[0 length(tides_out.lat)+1]);
    aa = axis;
    grid on;
    box on;
    set(gca,'xticklabel',[]);
    for kk = 1:length(tides_out.lat)
        txh = text(kk-.2,max(tides_out.(constname).amp(kk))+.12*diff(aa(3:4)),...
                   tides_out.station_name(kk),'fontsize',8);
        set(txh,'rotation',60);
    end
    clear kk
    % set the upper ylim to ~30% more than the max amplitude (to fit names)
    set(gca,'ylim',[0 ceil(10.*1.3.*(max(max(tides_out.(constname).amp))))./10]);
    ylabel('Amplitude [m]','fontsize',12);
    % add a title
    %title(['Tidal constituent:  ' upper(constname)],'fontsize',14);
    % make the title text 10% above the upper ylim
    ttx = text(aa(1)+0.5,1.1.*(ceil(10.*1.3.*(max(max(tides_out.(constname).amp))))./10),...
               ['Tidal constituent:  ' upper(constname)]);
    set(ttx,'fontweight','bold','fontsize',14);
    
    subplot_close(2,1,2); hold on;
    bph = bar(1:length(tides_out.lat),tides_out.(constname).pha);
    set(gca,'xlim',[0 length(tides_out.lat)+1]);
    aa = axis;
    axis([0 length(tides_out.lat)+1 0 360])
    grid on;
    box on;
    set(gca,'ytick',[0:60:360],'yticklabel',[0:60:360]);
    set(gca,'xtick',[1:length(tides_out.lat)]);
    set(gca,'xticklabel',round(10.*tides_out.lat)./10);
    ylabel('Phase [deg]','fontsize',12);
    xlabel('Latitude [deg]','fontsize',12);
    [LEGH,OBJH,OUTH,OUTM] = legend('obs','ROMS');
    set(LEGH,'box','off','orientation','horizontal','location','northeast');
end
clear ii constname


%% Load coastal-Salish tide tracks and compute along-track distance
load coastSalish_tide_track.mat % track_PS track_SG

% Calculate and store distances between individual adjacent lon,lat pairs (Puget Sound)
track_PSdist = [0]; % initialize with zero for first along-track distance
for ii = 1:length(track_PS)-1
    [dist, phase] = sw_dist([track_PS(ii,2) track_PS(ii+1,2)],...
        [track_PS(ii,1) track_PS(ii+1,1)],'km');
    track_PSdist = [track_PSdist; dist];
    clear dist phase
end
clear ii
track_PSdist = cumsum(track_PSdist); % along-track distance from first element

% Calculate and store distances between individual adjacent lon,lat pairs (Strait of Georgia)
track_SGdist = [0]; % initialize with zero for first along-track distance
for ii = 1:length(track_SG)-1
    [dist, phase] = sw_dist([track_SG(ii,2) track_SG(ii+1,2)],...
        [track_SG(ii) track_SG(ii+1)],'km');
    track_SGdist = [track_SGdist; dist];
    clear dist phase
end
clear ii
track_SGdist = cumsum(track_SGdist); % along-track distance from first element


%% Determine along-track distance to Neah Bay, WA (the entrance to JdF)
% So we can normalize track distance relative to JdF entrance
%
% Identify Neah Bay
tg_id = NaN.*ones([length(tides_out.lat) 1]);
for kk = 1:length(tides_out.lat)
    tg_id(kk) = strcmp(tides_out.station_name(kk),'Neah Bay');
end
clear kk
tg_id = find(tg_id == 1);
%
% Calculate distance from gauge to all points in the transect
transect2gauge = [];
for kk = 1:length(track_PS)
    [dist, phase] = sw_dist([track_PS(kk,2) tides_out.lat(tg_id)],...
        [track_PS(kk,1) tides_out.lon(tg_id)],'km');
    transect2gauge = [transect2gauge; dist];
    clear dist phase
end
clear kk
%
% Find the closest point on the transect to the current gauge
mindist_id = find(abs(transect2gauge) == min(abs(transect2gauge)));
clear transect2gauge
% Along-transect distance to Neah Bay, WA
atd_NBdist = track_PSdist(mindist_id);
clear mindist_id


%% Update track distances so they are relative to Neah Bay
track_PSdist = track_PSdist - atd_NBdist;
track_SGdist = track_SGdist - atd_NBdist;

% Calculate min and max values for xlim in along-transect distance plots
minx = 100.*floor(min(track_SGdist)./100);
maxx = 100.*ceil(max(track_SGdist)./100);


%% Set up desired tide gauges to plot
% Below we break them up into different sectors - coast, JdF, PS, SG
% (at minimum we need the PS, and SG as separate transects)
%
% Tide gague stations for the outer coast
tg_CST = {'Charleston'
         'South Beach'
         'Garibaldi'
         'Astoria'
         'Toke Point'
         'Westport'
         'La Push'
         'Neah Bay'};
% Tide gague stations for Juan de Fuca Strait
tg_JdF = {'Neah Bay'
         'Port Angeles'
         'Victoria'};
% Tide gauge stations to include for Puget Sound
tg_PS = {'Port Angeles'
         'Port Townsend'
         'Seattle'
         'Tacoma'};
% Tide gauge stations to include for Strait of Georgia
tg_SG = {'Victoria'
         'Friday Harbor'
         'Patricia Bay'
         'Cherry Point'
         'Vancouver'
         'Point Atkinson'
         'Campbell River'};


%% Gather tide information at the predefined coastal stations...
% Initialize storage
count = 0; % keeps track of number of active tide gauges (having data)
cstid_vec = []; % ID of tide gauges with data
atd_CSTlon = NaN.*ones(length(tg_CST),1); % lon in track closest to gauge
atd_CSTlat = NaN.*ones(length(tg_CST),1); % lat in track closest to gauge
atd_CSTdist = NaN.*ones(length(tg_CST),1); % approximate along-track distance of gauge

for jj = 1:length(tg_CST)
    
    % Identify if the current tide gauge (tg_CST) exists in the active list (TG)
    tg_id = NaN.*ones([length(tides_out.lat) 1]);
    for kk = 1:length(tides_out.lat)
        tg_id(kk) = strcmp(tides_out.station_name(kk),tg_CST(jj));
    end
    clear kk
    tg_id = find(tg_id == 1);
    
    % Proceed if the tide gauge in tg_PS exists in TG
    if ~isempty(tg_id)
        
        count = count+1; % update number of tide gauges with data
        cstid_vec = [cstid_vec; tg_id]; % store the tide gauge ID
        
        % Calculate distance from gauge to all points in the transect
        transect2gauge = [];
        for kk = 1:length(track_PS)
            [dist, phase] = sw_dist([track_PS(kk,2) tides_out.lat(tg_id)],...
                [track_PS(kk,1) tides_out.lon(tg_id)],'km');
            transect2gauge = [transect2gauge; dist];
            clear dist phase
        end
        clear kk
        
        % Find the closest point on the transect to the current gauge
        mindist_id = find(abs(transect2gauge) == min(abs(transect2gauge)));
        clear transect2gauge
        
        % Store the closest lon, lat, along-track distance on the transect
        atd_CSTlon(count) = track_PS(mindist_id,1);
        atd_CSTlat(count) = track_PS(mindist_id,2);
        atd_CSTdist(count) = track_PSdist(mindist_id);
        clear mindist_id
        
    end
    
    clear tg_id
end
clear jj count

% NaNs may exist because we pre-allocated the size too large
atd_CSTlon = atd_CSTlon(isfinite(atd_CSTlon));
atd_CSTlat = atd_CSTlat(isfinite(atd_CSTlat));
atd_CSTdist = atd_CSTdist(isfinite(atd_CSTdist));

     
%% Gather tide information at the predefined Juan de Fuca stations...
% Initialize storage
count = 0; % keeps track of number of active tide gauges (having data)
jdfid_vec = []; % ID of tide gauges with data
atd_JdFlon = NaN.*ones(length(tg_JdF),1); % lon in track closest to gauge
atd_JdFlat = NaN.*ones(length(tg_JdF),1); % lat in track closest to gauge
atd_JdFdist = NaN.*ones(length(tg_JdF),1); % approximate along-track distance of gauge

for jj = 1:length(tg_JdF)
    
    % Identify if the current tide gauge (tg_JdF) exists in the active list (TG)
    tg_id = NaN.*ones([length(tides_out.lat) 1]);
    for kk = 1:length(tides_out.lat)
        tg_id(kk) = strcmp(tides_out.station_name(kk),tg_JdF(jj));
    end
    clear kk
    tg_id = find(tg_id == 1);
    
    % Proceed if the tide gauge in tg_PS exists in TG
    if ~isempty(tg_id)
        
        count = count+1; % update number of tide gauges with data
        jdfid_vec = [jdfid_vec; tg_id]; % store the tide gauge ID
        
        % Calculate distance from gauge to all points in the transect
        transect2gauge = [];
        for kk = 1:length(track_PS)
            [dist, phase] = sw_dist([track_PS(kk,2) tides_out.lat(tg_id)],...
                [track_PS(kk,1) tides_out.lon(tg_id)],'km');
            transect2gauge = [transect2gauge; dist];
            clear dist phase
        end
        clear kk
        
        % Find the closest point on the transect to the current gauge
        mindist_id = find(abs(transect2gauge) == min(abs(transect2gauge)));
        clear transect2gauge
        
        % Store the closest lon, lat, along-track distance on the transect
        atd_JdFlon(count) = track_PS(mindist_id,1);
        atd_JdFlat(count) = track_PS(mindist_id,2);
        atd_JdFdist(count) = track_PSdist(mindist_id);
        clear mindist_id
        
    end
    
    clear tg_id
end
clear jj count

% NaNs may exist because we pre-allocated the size too large
atd_JdFlon = atd_JdFlon(isfinite(atd_JdFlon));
atd_JdFlat = atd_JdFlat(isfinite(atd_JdFlat));
atd_JdFdist = atd_JdFdist(isfinite(atd_JdFdist));

     
%% Gather tide information at the predefined Puget Sound stations...
% Initialize storage
count = 0; % keeps track of number of active tide gauges (having data)
psid_vec = []; % ID of tide gauges with data
atd_PSlon = NaN.*ones(length(tg_PS),1); % lon in track closest to gauge
atd_PSlat = NaN.*ones(length(tg_PS),1); % lat in track closest to gauge
atd_PSdist = NaN.*ones(length(tg_PS),1); % approximate along-track distance of gauge

for jj = 1:length(tg_PS)
    
    % Identify if the current tide gauge (tg_PS) exists in the active list (TG)
    tg_id = NaN.*ones([length(tides_out.lat) 1]);
    for kk = 1:length(tides_out.lat)
        tg_id(kk) = strcmp(tides_out.station_name(kk),tg_PS(jj));
    end
    clear kk
    tg_id = find(tg_id == 1);
    
    % Proceed if the tide gauge in tg_PS exists in TG
    if ~isempty(tg_id)
        
        count = count+1; % update number of tide gauges with data
        psid_vec = [psid_vec; tg_id]; % store the tide gauge ID
        
        % Calculate distance from gauge to all points in the transect
        transect2gauge = [];
        for kk = 1:length(track_PS)
            [dist, phase] = sw_dist([track_PS(kk,2) tides_out.lat(tg_id)],...
                [track_PS(kk,1) tides_out.lon(tg_id)],'km');
            transect2gauge = [transect2gauge; dist];
            clear dist phase
        end
        clear kk
        
        % Find the closest point on the transect to the current gauge
        mindist_id = find(abs(transect2gauge) == min(abs(transect2gauge)));
        clear transect2gauge
        
        % Store the closest lon, lat, along-track distance on the transect
        atd_PSlon(count) = track_PS(mindist_id,1);
        atd_PSlat(count) = track_PS(mindist_id,2);
        atd_PSdist(count) = track_PSdist(mindist_id);
        clear mindist_id
        
    end
    
    clear tg_id
end
clear jj count

% NaNs may exist because we pre-allocated the size too large
atd_PSlon = atd_PSlon(isfinite(atd_PSlon));
atd_PSlat = atd_PSlat(isfinite(atd_PSlat));
atd_PSdist = atd_PSdist(isfinite(atd_PSdist));


%% Gather tide information at the predefined Strait of Georgia stations...
% Initialize storage
count = 0; % keeps track of number of active tide gauges (having data)
sgid_vec = []; % ID of tide gauges with data
atd_SGlon = NaN.*ones(length(tg_SG),1); % lon in track closest to gauge
atd_SGlat = NaN.*ones(length(tg_SG),1); % lat in track closest to gauge
atd_SGdist = NaN.*ones(length(tg_SG),1); % approximate along-track distance of gauge

for jj = 1:length(tg_SG)
    
    % Identify if the current tide gauge (tg_SG) exists in the active list (TG)
    tg_id = NaN.*ones([length(tides_out.lat) 1]);
    for kk = 1:length(tides_out.lat)
        tg_id(kk) = strcmp(tides_out.station_name(kk),tg_SG(jj));
    end
    clear kk
    tg_id = find(tg_id == 1);
    
    % Proceed if the tide gauge in tg_SG exists in TG
    if ~isempty(tg_id)
        
        count = count+1; % update number of tide gauges with data
        sgid_vec = [sgid_vec; tg_id]; % store the tide gauge ID
        
        % Calculate distance from gauge to all points in the transect
        transect2gauge = [];
        for kk = 1:length(track_SG)
            [dist, phase] = sw_dist([track_SG(kk,2) tides_out.lat(tg_id)],...
                [track_SG(kk,1) tides_out.lon(tg_id)],'km');
            transect2gauge = [transect2gauge; dist];
            clear dist phase
        end
        clear kk
        
        % Find the closest point on the transect to the current gauge
        mindist_id = find(abs(transect2gauge) == min(abs(transect2gauge)));
        clear transect2gauge
        
        % Store the closest lon, lat, along-track distance on the transect
        atd_SGlon(count) = track_SG(mindist_id,1);
        atd_SGlat(count) = track_SG(mindist_id,2);
        atd_SGdist(count) = track_SGdist(mindist_id);
        clear mindist_id
        
    end
    
    clear tg_id
end
clear jj count

% NaNs may exist because we pre-allocated the size too large
atd_SGlon = atd_SGlon(isfinite(atd_SGlon));
atd_SGlat = atd_SGlat(isfinite(atd_SGlat));
atd_SGdist = atd_SGdist(isfinite(atd_SGdist));


%% Plot a map with the stations and tracks
figure; hold on;
set(gcf,'position',[1300 290 400 630]);
plot_WAcoast('all','color',0.6*[1 1 1]);
axis([-127 -122 43 50.5]);
Z_dar;

for ii = 1:length(tides_out.lat)
    plot(tides_out.lon(ii),tides_out.lat(ii),'ko','markerfacecolor','g','markersize',6);
    
    % label station names
    switch tides_out.station_name{ii}
        case {'Garibaldi','Toke Point','Westport','La Push','Port Townsend','Vancouver','Campbell River','Seattle','Tacoma'}
            text(tides_out.lon(ii)+0.2,tides_out.lat(ii)-0.15,tides_out.station_name(ii));
        case 'Neah Bay'
            text(tides_out.lon(ii)-1.2,tides_out.lat(ii)+0.1,tides_out.station_name(ii));
        case 'Point Atkinson'
            text(tides_out.lon(ii)+0.2,tides_out.lat(ii)+0.1,tides_out.station_name(ii));
        case 'Friday Harbor'
            text(tides_out.lon(ii)+0.2,tides_out.lat(ii)-0.05,tides_out.station_name(ii));
        case 'Victoria'
            text(tides_out.lon(ii)+0.1,tides_out.lat(ii)-0.07,tides_out.station_name(ii));
        case 'Astoria'
            text(tides_out.lon(ii),tides_out.lat(ii)-0.2,tides_out.station_name(ii));
        case 'Port Angeles'
            text(tides_out.lon(ii)-0.75,tides_out.lat(ii)+0.07,tides_out.station_name(ii));
        otherwise
            text(tides_out.lon(ii)+0.2,tides_out.lat(ii),tides_out.station_name(ii));
    end
end
xlabel('Longitude [deg]','fontsize',12);
ylabel('Latitude [deg]','fontsize',12);

% Add the tracks...
plot(track_PS(:,1),track_PS(:,2),'m-');
plot(track_SG(:,1),track_SG(:,2),'m-');
%print(gcf,'-dpng',[outdir 'TideGaugeSites_transects']);


%% Plot tidal amplitudes and phases as a function of distance from Neah Bay
for ii = 1:length(constit)
    constname = char(constit(ii)); % the current constituent
    
    figure; hold on;
    set(gcf,'position',[1000 300 900 450]);
    
    % --TIDAL AMPLITUDES--
    subplot_close(2,1,1); hold on;
    % Plot Puget Sound stations
    % observations
    pho_ps = plot(atd_PSdist,tides_out.(constname).amp(psid_vec,1),'go-',... 
        'markerfacecolor','g','markersize',7,'linewidth',2);
    % roms
    phm_ps = plot(atd_PSdist,tides_out.(constname).amp(psid_vec,2),'gs-',... 
        'markersize',7,'linewidth',1);
    % Plot Strait of Georgia stations
    % obs
    pho_sg = plot(atd_SGdist,tides_out.(constname).amp(sgid_vec,1),'bo-',... 
        'markerfacecolor','b','markersize',7,'linewidth',2);
    % roms
    phm_sg = plot(atd_SGdist,tides_out.(constname).amp(sgid_vec,2),'bs-',...
        'markersize',7,'linewidth',1);
    % Plot Coastal stations
    % obs
    pho_cst = plot(atd_CSTdist,tides_out.(constname).amp(cstid_vec,1),'ko-',... 
        'markerfacecolor','k','markersize',7,'linewidth',2);
    % roms
    phm_cst = plot(atd_CSTdist,tides_out.(constname).amp(cstid_vec,2),'ks-',... 
        'markersize',7,'linewidth',1);
    % Plot Juan de Fuca stations
    % obs
    pho_jdf = plot(atd_JdFdist,tides_out.(constname).amp(jdfid_vec,1),'ro-',... 
        'markerfacecolor','r','markersize',7,'linewidth',2);
    % roms
    phm_jdf = plot(atd_JdFdist,tides_out.(constname).amp(jdfid_vec,2),'rs-',... 
        'markersize',7,'linewidth',1);
    set(gca,'xlim',[minx maxx]);
    set(gca,'xtick',[minx:100:maxx]);
    set(gca,'xticklabel',[]);
    % set the upper ylim to ~20% more than the max amplitude (to fit names)
    set(gca,'ylim',[0 ceil(10.*1.2.*(max(max(tides_out.(constname).amp))))./10]);
    ylabel('Amplitude [m]','fontsize',12);
    aa = axis;
    grid on;
    box on;
    % add a title
    %title(['Tidal constituent:  ' upper(constname)],'fontsize',14);
    % make the title text 10% above the upper ylim
    ttx = text(aa(1)+0.5,1.1.*(ceil(10.*1.2.*(max(max(tides_out.(constname).amp))))./10),...
               ['Tidal constituent:  ' upper(constname)]);
    set(ttx,'fontweight','bold','fontsize',14);
    % Add a legend
    %[LEGH,OBJH,OUTH,OUTM] = legend([pho_cst, pho_jdf, pho_ps, pho_sg],...
    %    'Outer coast','Juan de Fuca','Puget Sound','Strait of Georgia');
    %set(LEGH,'box','off','orientation','horizontal','location','southeast');
    [LEGH,OBJH,OUTH,OUTM] = legend([pho_cst, phm_cst],'observations','model');
    set(LEGH,'box','off','orientation','vertical','location','southwest');
    % Add colored text
    tx_cst = text(aa(1)+0.2.*(aa(2)-aa(1)),aa(3)-0.1.*(aa(4)-aa(3)),'Outer coast');
    set(tx_cst,'fontweight','bold','color','k');
    tx_jdf = text(aa(1)+0.35.*(aa(2)-aa(1)),aa(3)-0.1.*(aa(4)-aa(3)),'Juan de Fuca');
    set(tx_jdf,'fontweight','bold','color','r');
    tx_ps = text(aa(1)+0.5.*(aa(2)-aa(1)),aa(3)-0.1.*(aa(4)-aa(3)),'Puget Sound');
    set(tx_ps,'fontweight','bold','color','g');
    tx_sg = text(aa(1)+0.65.*(aa(2)-aa(1)),aa(3)-0.1.*(aa(4)-aa(3)),'Strait of Georgia');
    set(tx_sg,'fontweight','bold','color','b');
    
    
    
    % Before moving to phase plots, identify southernmost coastal station
    cst_start_id = find(tides_out.lat(cstid_vec) == min(tides_out.lat(cstid_vec)));
    
    % Determine the initial phase to remove from the phase line plots (for plotting)
    clear initial_obs_pha initial_roms_pha initial_pha
    initial_obs_pha = tides_out.(constname).pha(cstid_vec(cst_start_id),1);
    initial_roms_pha = tides_out.(constname).pha(cstid_vec(cst_start_id),2);
    initial_pha = min([initial_obs_pha initial_roms_pha]);
    
    % Redefine phase as relative to southernmost site (for plotting)
    ps_obs_pha = (tides_out.(constname).pha(psid_vec,1)) - initial_pha;
    ps_roms_pha = (tides_out.(constname).pha(psid_vec,2)) - initial_pha;
    sg_obs_pha = (tides_out.(constname).pha(sgid_vec,1)) - initial_pha;
    sg_roms_pha = (tides_out.(constname).pha(sgid_vec,2)) - initial_pha;
    cst_obs_pha = (tides_out.(constname).pha(cstid_vec,1)) - initial_pha;
    cst_roms_pha = (tides_out.(constname).pha(cstid_vec,2)) - initial_pha;
    jdf_obs_pha = (tides_out.(constname).pha(jdfid_vec,1)) - initial_pha;
    jdf_roms_pha = (tides_out.(constname).pha(jdfid_vec,2)) - initial_pha;
    
    % Add 360 deg if necessary for wrapping of the phase vector
    for jj = 1:length(ps_obs_pha) % obs and roms should have same lengths
        if ps_obs_pha(jj) < 0
            ps_obs_pha(jj) = ps_obs_pha(jj) + 360;
        end
        if ps_roms_pha(jj) < 0
            ps_roms_pha(jj) = ps_roms_pha(jj) + 360;
        end
    end
    clear jj
    for jj = 1:length(sg_obs_pha) % obs and roms should have same lengths
        if sg_obs_pha(jj) < 0
            sg_obs_pha(jj) = sg_obs_pha(jj) + 360;
        end
        if sg_roms_pha(jj) < 0
            sg_roms_pha(jj) = sg_roms_pha(jj) + 360;
        end
    end
    clear jj
    for jj = 1:length(cst_obs_pha) % obs and roms should have same lengths
        if cst_obs_pha(jj) < 0
            cst_obs_pha(jj) = cst_obs_pha(jj) + 360;
        end
        if cst_roms_pha(jj) < 0
            cst_roms_pha(jj) = cst_roms_pha(jj) + 360;
        end
    end
    clear jj
    for jj = 1:length(jdf_obs_pha) % obs and roms should have same lengths
        if jdf_obs_pha(jj) < 0
            jdf_obs_pha(jj) = jdf_obs_pha(jj) + 360;
        end
        if jdf_roms_pha(jj) < 0
            jdf_roms_pha(jj) = jdf_roms_pha(jj) + 360;
        end
    end
    clear jj
    
    
    
    % --TIDAL PHASES--
    subplot_close(2,1,2); hold on;
    % Plot Puget Sound stations
    % obs
    plot(atd_PSdist,ps_obs_pha,'go-','markerfacecolor','g','markersize',7,'linewidth',2);
    % roms
    plot(atd_PSdist,ps_roms_pha,'gs-','markersize',7,'linewidth',1);
    % Plot Strait of Georgia stations
    % obs
    plot(atd_SGdist,sg_obs_pha,'bo-','markerfacecolor','b','markersize',7,'linewidth',2);    
    % roms
    plot(atd_SGdist,sg_roms_pha,'bs-','markersize',7,'linewidth',1);
    % Plot Coastal stations
    % obs
    plot(atd_CSTdist,cst_obs_pha,'ko-','markerfacecolor','k','markersize',7,'linewidth',2);
    % roms
    plot(atd_CSTdist,cst_roms_pha,'ks-','markersize',7,'linewidth',1);
    % Plot Juan de Fuca stations
    % obs
    plot(atd_JdFdist,jdf_obs_pha,'ro-','markerfacecolor','r','markersize',7,'linewidth',2);
    % roms
    plot(atd_JdFdist,jdf_roms_pha,'rs-','markersize',7,'linewidth',1);
    
%  Commented out section is the original much simpler version where we 
%  simply plot phases as they are instead of relative to southernmost site...
%     % Plot Puget Sound stations
%     % obs
%     plot(atd_PSdist,(tides_out.(constname).pha(psid_vec,1)),'go-',... 
%         'markerfacecolor','g','markersize',7,'linewidth',2);
%     % roms
%     plot(atd_PSdist,(tides_out.(constname).pha(psid_vec,2)),'gs-',... 
%         'markersize',7,'linewidth',1);
%     % Plot Strait of Georgia stations
%     % obs
%     plot(atd_SGdist,(tides_out.(constname).pha(sgid_vec,1)),'bo-',... 
%         'markerfacecolor','b','markersize',7,'linewidth',2);
%     % roms
%     plot(atd_SGdist,(tides_out.(constname).pha(sgid_vec,2)),'bs-',... 
%         'markersize',7,'linewidth',1);
%     % Plot Coastal stations
%     % obs
%     plot(atd_CSTdist,(tides_out.(constname).pha(cstid_vec,1)),'ko-',... 
%         'markerfacecolor','k','markersize',7,'linewidth',2);
%     % roms
%     plot(atd_CSTdist,(tides_out.(constname).pha(cstid_vec,2)),'ks-',... 
%         'markersize',7,'linewidth',1);
%     % Plot Juan de Fuca stations
%     % obs
%     plot(atd_JdFdist,(tides_out.(constname).pha(jdfid_vec,1)),'ro-',... 
%         'markerfacecolor','r','markersize',7,'linewidth',2);
%     % roms
%     plot(atd_JdFdist,(tides_out.(constname).pha(jdfid_vec,2)),'rs-',... 
%         'markersize',7,'linewidth',1);
    
    set(gca,'xlim',[minx maxx]);
    set(gca,'xtick',[minx:100:maxx]);
    set(gca,'xticklabel',[minx:100:maxx]);
    set(gca,'ylim',[0 360]);
    set(gca,'ytick',[0:60:360],'yticklabel',[0:60:360]);
    box on;
    grid on;
    ylabel('*Phase [deg]','fontsize',12);
    xlabel('Inland distance from Neah Bay, WA [km]','fontsize',12);
    % disclaimer!
    tx_dc = text(minx+50,310,'*Phase relative to southernmost coastal site');
    set(tx_dc,'fontsize',fs1-4);
    
end
clear ii constname
